﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary;
using SmartMathLibrary.LinearEquationSolvers.Preconditioners;

namespace SmartMathLibrary.LinearEquationSolvers
{
    /// <summary>
    /// This class provides the solving of a linear system of equation by using the RichardsonIteration method.
    /// The algorithm, which will be used partly based on the IML++ library.
    /// </summary>
    [Serializable]
    public class RichardsonIterationSystemSolver : AbstractNonStationaryIterativeMethod
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="RichardsonIterationSystemSolver"/> class.
        /// </summary>
        /// <param name="coefficientMatrix">The coefficient matrix of the equation system to solve.</param>
        /// <param name="resultVector">The result vector of the equation system to solve.</param>
        /// <param name="preconditioner">The preconditioner which should be used for the iterative method.</param>
        public RichardsonIterationSystemSolver(Matrix coefficientMatrix, GeneralVector resultVector,
                                               IPreconditioner preconditioner)
            : base(coefficientMatrix, resultVector, preconditioner)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RichardsonIterationSystemSolver"/> class.
        /// </summary>
        /// <param name="equationSystem">The equation system to solve.</param>
        /// <param name="preconditioner">The preconditioner which should be used for the iterative method.</param>
        public RichardsonIterationSystemSolver(LinearEquationSystem equationSystem, IPreconditioner preconditioner)
            : base(equationSystem, preconditioner)
        {
        }

        /// <summary>
        /// Solves the linear system of equation.
        /// </summary>
        /// <param name="startValues">The start values for x.</param>
        /// <param name="iterations">The maximum numbers of iterations to solve the system of equation.</param>
        /// <param name="precision">The precision of the result.</param>
        /// <returns>
        /// The vector which solves the linear system of equation Ax = b.
        /// </returns>
        public override GeneralVector SolveLinearSystem(GeneralVector startValues, int iterations, double precision)
        {
            //The algorithm partly based on the IML++ library of the National Institute of Standards and Technology.
            double residual;
            double lengthB = this.ResultVector.Length;
            GeneralVector r = this.ResultVector - this.CoefficientMatrix * startValues;
            Matrix preconditionMatrix = this.Preconditioner.ComputePreconditionMatrix();

            if (lengthB == 0.0)
            {
                lengthB = 1;
            }

            if ((residual = r.Length / lengthB) <= precision)
            {
                this.relativeError = residual;
                this.neededIterations = 0;
                this.precisionError = false;

                return startValues;
            }

            for (int i = 0; i < iterations; i++)
            {
                GeneralVector z = preconditionMatrix * r;

                startValues += z;
                r = this.ResultVector - this.CoefficientMatrix * startValues;

                if ((residual = r.Length / lengthB) <= precision)
                {
                    this.relativeError = residual;
                    this.neededIterations = i + 1;
                    this.precisionError = false;

                    return startValues;
                }
            }

            this.relativeError = residual;
            this.neededIterations = iterations;
            this.precisionError = true;

            return startValues;
        }
    }
}